Scaling in the Rubinstein-Duke Model for Reptation 
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We consider an arbitrarily charged polymer driven by a weak field through a gel according to the 
rules of the Rubinstein-Duke model. The probability distribution in the stationary state is related 
to that of the model in which only the head is charged. Thereby drift velocity, diffusion constant and 
orientation of any charged polymers are expressed in terms of those of the central model. Mapping 
the problem on a random walk of a tagged particle along a one-dimenional chain, leads to a unified 
scaling expression for the local orientation. It provides also an elucidation of the role of corrections 
to scaling. 
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The basic ingredient of the physics of polymers is scal- 
ing, i.e. the analysis of the properties as function of 
the length N of the polymer chain l|. While most of 
the leading scaling behavior is understood, the compari- 
son between theory and experiment is hampered by large 
corrections to scaling. A typical example is the renewal 
time t for polymers dissolved in a gel. Theory predicts 
behavior r ~ A 3 , while experiments seem to converge 
on r ~ N 3A Q. The discrepancy has been blamed on 
large corrections to scaling, blurring the true asymptotic 
behavior 0, [j| . This suggestion has been given a firm 
basis by the analysis of Carlon et. al. 4] of the re- 
newal time in the Rubinstein-Duke (RD) model for rep- 
tation. Polymers up to thousands of base pairs may still 
have large corrections to the asymptotic scaling behavior. 
Since they find that deviations from the leading scaling 
behavior decay as 1/y/N, their work indicates that the 
correction to scaling exponent is -1/2. 

Sofar the scaling analysis has been restricted to global 
properties, such as renewal time and diffusion coefficient. 
Equally interesting and more informative, are scaling 
properties related to the position in the chain. As ex- 
ample consider a neutral polymer dissolved in a gel with 
a magnetic beat attached to the head of the polymer. 
The polymer can be pulled through the gel by a mag- 
netic field. It is the analogon of the more common elec- 
trophoresis (EP) and has been named magnetophoresis 
(MP) [jj. The field orients the polymer links, at the head 
stronger than at the tail. On top of this overal effect is 
a subtile scaling behavior near head and tail, as we will 
show. Similar effects occur in electrophoresis, which is 
a central tool of DNA fingerprinting. DNA, being an 
acid, acquires as local charge, when dissolved in a gel 
like agarose. These local charges are evenly spread over 
the chain and an electric field pulls equally strong at all 
elements of the polymer. The local orientation of this 
chain has an even more interesting scaling behavior. 

The motion of a polymer is far too complicated to be 



taken into full detailed account. Therefore lattice models 
have been designed for reptation, in which the polymer 
is viewed as a chain of hopping reptons [6|, confined to 
a tube of pores in the gel. Among them the RD model 
stands out by the simplicity of the motion rules 0, 0, • 
Although the model is a rather crude simplication of the 
reality, it captures a number of basic features of reptation 
. This paper is concerned with scaling effects involving 
the position along the chain, as they emerge in the RD 
model. 

An attractive aspect of the RD model is the fact that 
its hopping operator is a one-dimensional spin operator, 
albeit a non-hermitian hamiltonian. This allows to ap- 
ply the Density Matrix Renormalization Group (DMRG) 
technique to study the stationary state of the Master 
Equation 0, The DMRG method gives precision 

data for the whole scala of chains up to lengths of 150 
reptons, which make it ideally suited for finite size analy- 
sis. This has an advantage over direct simulations of the 
model, which are slow, due to the long renewal time, and 
limited by statistical errors. 

In this paper we linearize the Master Equation for the 
probability distribution with respect to the driving field 
and study the solution for an arbitrarily charged poly- 
mer. First we relate the probability distribution of the 
general case to that for the MP model. This enables us to 
link the properties of the original RD model, with equal 
charges on all reptons, to those of the MP model. We 
map the MP model on the problem of a random walker 
on a one-dimensional chain. Going back and forth from 
the EP and MP variant, we derive scaling properties for 
the various regions of the local orientation. 

The RD model views the polymer as a a one- 
dimensional string of reptons connected by links. The 
reptons are located in the cells of a rf-dimensional hyper- 
cubic lattice with the field along the body diagonal 
We number the reptons from 1 (tail) to N (head). The 
reptons hop independently from each other, under the 
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constraint that the links always connect two reptons in 
adjacent cells, or two reptons in the same cell. Only the 
projection of a link on the field direction is important for 
the probability distribution. A configuration of the poly- 
mer can be represented by a vector y = (yi, ■ ■ ■ , yjv-i), 
of link variables, with yj — 0, ±1, measuring the distance 
between the reptons j and j + 1 along the field direction 
0. 

The Master Equation for the stationary probability 
distribution P(y) has the form 

£ [W(y\y')P(y') - W(y'\y)P(y)} = MP(y) = 0. (1) 
y 

W(y|y') is the transition rate from configuration y' to 
y. A move of repton i in the direction of the field has a 
bias Bi = exp(<^e/2), where e measures the strength of 
the field and qi the charge of the repton. Reptons mov- 
ing opposite to the field are biased with B~ x . Although 
the influence of the embedding dimension has interesting 
aspects 11], we confine ourselves in this paper to d = 1, 
for which the most extensive DMRG results are obtained, 
yielding the most accurate scaling data. 

We consider the field strength e as a small parameter 
(eN < 1); it is experimentally the most relevant regime. 
For e = the stochastic matrix is symmetric and has 
the solution P°(y) — i^^^ 1 ^, all S^ -1 configurations 
being equally probable. Expanding the Master Equation 
in powers of e 

M =M° + eM 1 + -- - , P(y) = P°(y) +eP 1 (y) + • • • , 

(2) 

leads to the lowest order equation Ai°P° = and to the 
first order equation 

M°P 1 +M X P° = 0. (3) 

This paper is concerned with its solution. 

The remarkable point is that the dependence on the 
charge distribution can be made explicit. 

N-l / j \ / N \ 

= - E E * %p0(y) + E 9i) pl Mp(yy 
j=i \i=i / \i=\ ) 

(4) 

Here P\ 1P (y) is the first order distribution of the MP 
model with a unit charge on the last repton only. The 
proof follows by substitution in ©. So all cases are re- 
duced to this central model. 

A direct consequence of this relation is that the drift 
velocity v of an arbitrary qi is related to that of the MP 
model as 

v = ^E^ vmp ( 5 ) 

Thus the drift velocity depends only on the total charge. 
This was anticipated in 12]. E.g. the drift velocity of 
the EP model is N times that of the MP model. 



A nice numerical illustration of relation (0J is given by 
the calculation of the local orientation (yj). This is a 
first order effect and (0} gives a relation e.g. between the 
EP model (all q t = 1) and the MP model (only q N = 1 ) . 
Both models have been treated with DMRG jllll l ll ^ . 
In Fig. ^we show the values for the MP model and the 
equivalent combination for the EP model. The corre- 
spondence is perfect. It shows that for both calculations 
e is small enough to guarantee that the weak field limit 
applies. It also is a proof of the accuracy of the DMRG 
method, since the calculations were performed prior to 
the derivation of relation 0}. Although the MP orien- 




FIG. 1: Comparison of the MP orientation with a translated 
EP values for the same value iV = 63 and d — 1 

tation is our basic ingredient, plotting the EP orienta- 
tions more clearly shows the N dependence, see Fig. El 
To explain the complicated scaling behavior of the these 
curves is the main aim of this paper. We already note 
that the curves (almost) pass through the same point at 
the head and tail, since there is a well known relation 
[fH IT3 ] between the orientation of the first link and the 
drift velocity. 

(yx) = -2e/3 + v EP = -2e/3 + e/3N. (6) 

In first order, the field dependence only combines with 
the sign of one link at the time 

JV-l 

Pmp(y)= E^W- ( ? ) 

Here x is the vector (ii, • • • ,xn-i) with Xj — The 
remaining factor Pj (x) depends on the position j of this 
link but not on the signs of the other links. To distin- 
guish link j from the others we give it a tag. The nota- 
tion Pj (x) anticipates that it can be seen as the (unnor- 
malized) probability that the chain is in configuration x 
and the tagged link at j. Relation Ql reduces the rel- 
evant configuration space from the original 3 N ^ X points 
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FIG. 2: Orientations for EP chains of length N = 31 to 151 
at e = 0.00005 and d = 1 



to (N — 1)2 N ^ 1 points (j,x). We refer to x% = 1 as a 
"particle" and to xi = as a "hole" . The marked link j 
is a tagged particle. 

The probability interpretation for Pj(x) is unique for 
the MP model. As J1J shows, all P 1 can be cast in the 
form Q, but the corresponding Pj's do not need to be 
positive. The equation for P,(x) follows from ©, sum- 
ming over all y which lead to the same configuration 



(J, x). The known term in J2J becomes with L = ^ 



j *3> 



M l MP P°{y) = yfr,N-i P°(x), P°(x) = 2 i /3 JV - 1 , 

(8) 

as only the last repton is effected by the field. We can 
eliminate this term from J2J using a symmetry of the EP 
model. When all reptons are pulled evenly, the transfor- 
mation 



y = (2/1, • • • , Vn-x) <-> y T = (-2/JV- 



-I/i), (9) 



leaves the probability invariant, as the numbering from 
head to tail is equivalent to our numbering from tail to 
head. This symmetry is clearly reflected in the curves 
of Fig. With relation @ we can transfer the EP 
symmetry to the MP model. It yields 



P J (x)+P A r_,(x T )=P°(x). 



(10) 



The transposed vector is x T = (x^r-i, ■ ■ ■ , aci). We use 
(|1(JII for j = N — 1, in order to eliminate the known term 
in @. Then equation J2J) for P, (x) becomes homoge- 
neous again and therefore a stationary Master Equation 
for Pj(x), which proves that Pj(x) can be interpreted 
as a probability. This new Master Equation is not very 
different from the original one QJ. The substitution has 
added an extra process: the tagged particle may jump 
from the tail to the head, thereby changing the configu- 
ration (1, x) to (N — 1, x T ). In addition we have the rule 
that the tagged particle may not leave or enter the chain, 
as a consequence of the fact that the tagged particle al- 
ways has to be present. 



This interpretation describes P/(x) as a random walk 
of a tagged particle at j in the sea of other particles 
and holes in configuration x. The possible moves are 
interchanges of particles and holes. Particles and holes 
may leave and enter the chain, but not the tag, which 
can jump from tail to head (and not reverse) as de- 
scribed above. The probability pj that the tag is at j, is 
^ x Pj(x). This is the quantity to be discussed, since it 
is equivalent to the local orientation: (yj) = 2pj/3, the 
2/3 coming from the fraction of particles in the chain. 

An equation for pj follows from the observation that 
the tag can only jump to neighboring positions j ± 1, 
leading to the effective equation 

W+(j - + W-(j + l) Pj+ i = [W+(j) + W-(j)]Pj, 

(11) 

with appropriate changes at the end links. W+(j) is the 
probability that, if the tag is at j, a hole is at j + 1 such 
that the tag can jump to the right. Similarly W~ (j) is the 
conditional probability that the tag at j is neighbored by 
a hole at j — 1. Fortunately the DMRG calculations not 
only provide information on pj , but also on the transition 
rates W±(j), since W±(J)pj is the joint probability on a 
tag-hole pair and these correlation functions have been 
calculated 0. If the holes were randomly distributed, 
one would expect 1 /3 for these rates, a value to which the 
rates converge in the bulk of the chain. The deviations 
from 1/3 are the largest at the tail of the chain (not at the 
head!). After a short initial regime they start to fall off 
smoothly as l/j. In the smooth regime one may replace 
the difference equation JTTJ, by a Fokker-Planck equa- 
tion, keeping only the first and second derivative with 
respect to j. 

The essential point in the scaling analysis is that the 
difference of the M^'s obeys a simple scaling relation 

[W+(j) - W-U)] - g(x)/j, x = j/VN, (12) 

which is demonstrated in Fig. [21 where we plot the dif- 
ference against the scaling variable x — j/yN. This 
difference enters in the Fokker-Planck equation as the 
systematic force. It tends to push the tag away from 
the tail. For N — > 00, the part at the head, where the 
various values of N fan out, shifts to larger and larger 
values, while the function g(x) approaches the limit 1/3. 
The small region at the tail side, where the scaling does 
not hold, shrinks with increasing N in this plot. In this 
region the transition rates change appreciably with each 
link, and one has to use the discrete equations ljll|l , start- 
ing from a value p\ . 

For values j ~ V~N the Fokker-Planck equation reads, 
writing pj — p(j/vN) 



d_ 

dx 



9{x) 



p{x) 



1 d 2 p(x) 
3 dx 2 



(13) 



In the right hand side we have replaced the sum of the 
W's by 2/3, since the correction is an order ViV smaller. 
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FIG. 3: The scaling form for the difference 1 12t 



not univocally find an exponent by fitting the initial ori- 
entation of a long chain. Here we see that the power 3 is 
fixed by the asymptotic properties of the scaling function 
g(x), for small and large x, following from vmp and the 
slope in the middle of the chain. We also note that this 
power does not set in inmediately at the tail, but devel- 
ops in the early regime of the variable x — j/yN. Even 
for our longest chains N = 151, the window where the 
power applies is rather small. 



In summary, the regions of order V N at the ends of 
the chain which differ from the bulk, originate from the 
algebraic decay (as towards the bulk of the corre- 
lations between two successive links. This explains why 
finite size corrections decay with an exponent —1/2. 



The dominant solution of Q13(l has the form 

p(x) =p(x )e R ( x \ R(x) = 3 f dx'g(x')/x'. (14) 



xq is the point where we match the solution with the 
outcome of the intial regime. 

The solution <|14ft bridges the initial and asymptotic 
behavior. We first inspect what (|14|> gives for the middle 
of the chain x m = >/~N jl. From the asymptotics of g{x) 
for small and large x one finds 



p(x m )^ P (x Q )N 3 ^+9(oo)]_ 



(15) 



Now p(x m ) = 1/2 by symmetry and p(xq) inherits its 
magnitude frompi, which equals 3ump/2 and is therefore 
order N~ 2 (see © and [3[^). Thus we find for the 
exponent 



3[ 3 (0) + s(oc)] 



(16) 



The value of g(oo) follows from the behavior of R(x) 
around x m . Expansion of g(x) in powers of x^ 1 gives 



n) 



3ff(oo) 



1 + 3 5 _i 



N 



N 



J 



(17) 



Thus 3(00 ) = 1/3, otherwise the behavior would not be 
linear in j/N. With JTJJ this implies 3(0) = 1. The 
1 / \/~N terms lead to the slope of the EP profile in the 
middle. This slope must grow as y/~N and a numerical fit 
gives g_i = 0.5. The small x behavior of R(x) yields 



p(x) ~ p(x ) (x/xqT 



(18) 



This power is higher than the estimate of Barkema and 
Newman[l6j|. who give 2.7 for the exponent. We could 
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